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Abstract. In this paper we present a new approach to computing 
homology (with field coefficients) and persistent homology. We 
use concepts from discrete Morse theory, to provide an algorithm 
which can be expressed solely in terms of simple graph theoretical 
operations. We use iterated Morse decomposition, which allows 
us to sidetrack many problems related to the standard discrete 
Morse theory. In particular, this approach is provably correct in 
any dimension. 



I. Preview. 

In this section we outline the purpose of the paper, assuming reader's 
familiarity with some concepts from computational topology. All the 
concepts will be carefully explained later. In this paper we introduce a 
new method to compute homology and persistent homology over field 
coefficients. The method is based on discrete Morse theory and is 
designed to be graph-theoretic. 

We present a brief, intuitive illustration of our method. As an ex- 
ample, let us consider a triangulation of a Dunce hat presented in 
Figure [Tja. We want to remind that this space has trivial homology 
but nontrivial topology. We will use discrete Morse theory to simplify 
the space, while preserving the homology. First, let us build a discrete 
Morse matching on this triangulation. It is well known that a Dunce hat 
has no perfect^ Morse complex [I]. Therefore, for any Morse matching 
we obtain some critical cells not corresponding to homology generators. 
The matching presented in Figure [T]b is optimal i.e. there few 
critical cells as possible. Now the Morse boundary is computed using 
the V-paths marked in Figure [TJc. The resulting Morse complex (with 
Z2 coefficients) is shown in Figure [TJd. Normally, in order to compute 
homology of a chain complex, boundary matrix is produced and Smith 



*A Morse complex is perfect if each critical cell corresponds to exactly one ho- 
mology generator. 
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Normal Form diagonalization is performed. However we would like to 
introduce an alternative approach. 

Using Kozlov's version of discrete Morse theory, we can iterate the 
Morse complex construction. In other words, we build a Morse com- 
plex of a Morse complex etc. So, here we compute a Morse matching 
of a chain complex presented in Figure [TJd. The only matched pair 
is marked with an arrow. Once the (iterated) Morse complex is com- 
puted, we are left with a single O-dimensional cell. It cannot be paired 
anymore and corresponds to the only homological feature: the con- 
nected component. This way we have computed the homology of the 
Dunce hat. 

Later the presented technique will be referred to as iterated Morse 
complex construction or iterated Morse decomposition. In this paper 
we will show that with the presented technique one can always acquire 
homology with field coefficients. We will also generalize it to the set- 
ting of persistent homology. As a result, we introduce a novel way to 
compute homology and persistence over a field. 



While persistent homology can potentially be applied to a plethora of 
different practical problems, ranging from sensor networks [35] to root 
architecture analysis [10], performance tends to be a problem. In par- 
ticular, there is growing interest in analysis of high- dimensional topo- 
logical features (going beyond connected components and 1-cycles). In 
low dimensions there exist efficient algorithms, but higher dimensional 
cases are still challenging. 

Such applications include analysis of complex networks such as social- 
network and biological networks such as gene-regulatory networks. Com- 
puting homology or persistent homology of maps between spaces leads 
to high dimensional datasets as the studied space is the Cartesian prod- 
uct of the source and target space [19]. 




The only class of algorithms to compute persistent homology in the 
general case is based on matrix reductions. Such an approach was 
recently shown to expose roughly quadratic computational complexity, 
for data coming from certain practical applications |38j. Also, it is 
not very suitable for distributed computing, which is a necessity as the 
datasets grow larger. 

Our aim is to propose an algorithmic framework which would scale 
reasonably well as the size of data (and its dimension) grows. In 
lower dimensions several algorithms were proposed and their efficiency 
stemmed from using fast techniques from graph-theory. Prompted by 
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Figure 1 . Iterated version of Morse complex construc- 
tion on a Dunce hat. On the top, the and 1 dimensional 
critical cells are marked with bold, the middle gray trian- 
gle is the unique critical 2-cell. In the bottom left picture 
the V-paths used to compute the Morse complex after 
the first iteration of the construction. On the bottom 
right, the second (and final) iteration of Morse complex 
construction. The remaining vertex corresponds to the 
unique homology generator in dimension 0. 



this observation, we wanted to propose a similar approach which would 
work for any dimension. 

However, directly extending the existing approaches to higher di- 
mensions is (as we believe) impossible. Combinatorial techniques, for 
instance described in [3, HH] , crucially depend on discrete Morse the- 
ory. In particular, a perfect Morse complex is (implicitly) constructed 
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in most cases, which allows to read homology information right away. 
Computing such perfect Morse complexes is hard (or even impossible) 
in higher dimensions. 

In case of field coefficients, which are important for practical applica- 
tions, the situation is more tractable. We use the properties of iterated 
Morse complexes [38], which always exist and are easy to compute. 
Additionally, we build on top of a recent theoretical framework by 
Mischaikow and Nanda [29J , which extends Morse theory to nitrations 
of spaces. 

The following paper describes the theory and algorithms for comput- 
ing homology and persistent homology using iterated Morse decompo- 
sition. We prove that the algorithm is correct for chain-complexes of 
any dimension. This includes commonly used simplicial and cubical 
complexes. Further, we show that, just as we intended, the algorithms 
can be entirely expressed in terms of basic graph-theoretical techniques. 
It promises that in terms of implementations, using efficient graph li- 
braries [27J will result in a scalable solution. 

The paper is structured as follows. Section [3] is a survey of exist- 
ing work in the topic. In Section [4] an introduction to homology and 
persistent homology is given, together with the necessary algorithmic 
background. In Section [5] a Discrete Morse Theory is highlighted. In 
Section [6] the concept of iterated Morse complex, crucial for this paper, 
is introduced. It is also explained there how this concept is used to com- 
pute homology. In Section [JJ these results are extended to simplification 
of filtered complexes, which is a preprocessing step for computing per- 
sistence. Later in Section [8] it is explained how to compute persistence 
using solely iterated Morse complex. Finally in Section [9] conclusions 
are drawn. 

3. Previous work. 

Computations of homology and persistent homology is a well estab- 
lished area of research with a rich history. In this section we want 
to summarize the main contributions and historical landmarks in this 
subject. 

The classical way of computing homology is by using Smith Normal 
Form (SNF) of a boundary operator matrix, see [30]. The classical 
algorithm has hyper-cubical complexity (in case of integer coefficients), 
see [36]. It is however possible to perform SNF in cubic time when field 
coefficients are considered. 

In the nineties Delfinado and Edelsbrunner provided an incremental 
algorithm for Betti numbers computation [TJ. This algorithm exhibits 
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linear time complexity and works for sub-triangulations of 3 dimen- 
sional sphere. There are also algorithms to compute homology with 
chain contractions and so called AT-models |15j . 

In 2003, after a few earlier iterations ( [39 1 IT4" 1 133] ) . persistent homol- 
ogy was introduced in its contemporary form [TT] . In the same paper, a 
matrix-reduction algorithm to compute persistence (Algorithm [T| was 
given. This is what is considered the standard algorithm to compute 
persistence. For a comprehensive presentation the reader should con- 
sult JTU]. Algorithmic results are further discussed in Section |4j 

Discrete Morse theory (DMT) was introduced by Robin Forman [13j . 
Later a more general, algebraic version was developed [25]. The com- 
prehensive presentation of algebraic discrete Morse theory can be found 
in [25]. The idea of using DMT for homology computations has been 
introduced by Lewiner [26] . The complexity aspects of DMT has been 
discussed in [22] ■ The notions of F— perfect and F— optimal Morse 
complexes are discussed in [TJ. A simplification algorithm for a Morse 
complexes on 2-manifolds has been presented in [2] . A divide and con- 
quer algorithm to compute Morse complexes has been presented in [17] . 

Recently Robins, Wood and Sheppard have provided a practical link 
between discrete Morse theory and persistence [31]. In this paper they 
introduce an optimal simplification scheme for persistence in case 3- 
dimensional complexes. The optimality of the presented result is re- 
stricted to 3 dimensions due to some deep results from simple homotopy 
theory. An extension of Morse theory suitable for simplifying filtered 
chain complexes in any dimension was later provided by Mischaikow 
and Nanda [29J. In short, by performing Morse matchings indepen- 
dently for each filtration level, a simplified filtered Morse complex is 
obtained. 

The idea of iterating Morse complex construction has already been 
used in [381 EE] as a tool to decrease the size of complexes before stan- 
dard algebraic computations. 

There are many software libraries to compute homology and per- 
sistent homology. For a homology software the reader should con- 
sult [31 El 121] ■ For persistent homology [231 El 132] are recommended. 

4. Background 

4.1. Complexes. We assume that the input data is represented as a 
chain complex with field coefficients. In the most typical case, this 
chain complex comes from a CW-decomposition of a given space which 
is a decomposition of a space into cells of different dimensions. In 
practice, simplicial and cubical complexes are used. For simplicity, we 
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will use Z 2 coefficients throughout the paper, as this is the standard 
setting for persistence. However, we want to remark that the presented 
algorithms work for any field coefficients. 

4.2. Boundary maps. Let us fix a complex /C. Cells of K, have 
different dimensions and are connected by boundary relations. If a 
(p — l)-cell a has non-zero boundary coefficient with a p-cell B, we say 
a is a proper face of B, and B is a proper coface of a. (Notation: cap- 
ital letters denote higher dimensional cell where a cell and its face is 
considered). Let a p- chain be a formal sum of p-cells with the Z 2 coeffi- 
cients. The boundary operator d p maps p-chains into p — 1-dimensional 
boundary chains. The chain of (co-)faces is called a (co-)boundary. We 
can extend the boundary operator linearly to p-chains. For any p-chain 
c = a i c ii we have d p c = a,id p Ci. It is assumed that the boundary 
of a boundary is zero, or formally: d p d p+ \ = 0. The p-chains, together 
with addition modulo 2, form a group of p- chains, denoted by C p . 

The boundary operator d p can be written as a binary matrix (also de- 
noted d p ), whose columns represent the boundaries and rows represent 
coboundaries of cells. 

4.3. Standard homology. Intuitively, homology can be used to cap- 
ture holes of complex /C. In 3-dimensional case holes are: connected 
components, tunnels, and voids. To define it formally, let us first in- 
troduce the group of p-cycles, Z p (fC) = kerd p and its subgroup: the 
group of p-boundaries, B p {fC) = imd p+ i. The p-th homology group is 
the quotient H p = Z p {fC) / B p {fC). The p-th. Betti number, denoted by 
P , is the rank of this group and counts the number of p-dimensional 
holes. 

4.4. Filtrations and persistence. For a given complex /C, a filtration 
is defined as a nested sequence of its subcomplexes: = /C_i C /Co C 
K\ C . . . C JC n — K [TO]- In case of persistence, filtrations are often 
generated by a filtering function, g : K — > Z defined on the input 
complex. We require that g(a) < g(B) whenever a is a face of B. This 
property guarantees that the sub-level sets JC t = 5 ,_1 (— oo,t] are sub- 
complexes of K, for each value of t 6 Z. The inclusions from /Q to )Cj, 
for i < j induce homomorphisms, / lJ : H(JCi) —> H{fCj). Complex K 
with filtration will be refered to as filtered complex. 

Given a complex K. and a filtering function g : K, — > Z, persistent 
homology studies homological changes of the sub-level complexes, K t = 
g~ l {— oo,t\. Persistent homology captures the birth and death times 
of homology classes of the sub-level complexes, as t grows from — oo to 
+oo. By birth, we mean that a homology feature is created; by death, 
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we mean it either becomes trivial or becomes identical to some other 
class born earlier. The persistence, or lifetime of a class, is the difference 
between the death and birth times. Often a multiset of persistence 
intervals is used to represent persistence. An interval encodes a lifetime 
of a homology class of a given dimension. We say that two spaces have 
the same persistence, if their corresponding persistence intervals are 
the same. 

The formal definition is as follows (after [10j): The p-th persistent 
homology groups of filtered complex K. are the images of the homo- 
morphisms induced by inclusion, H h ^{K) = imf l 'K For a standard 
definition of persistence diagram and persistence intervals the reader 
should consult [10J. 

We want to remind a theorem saying when persistence of two filtered 
complexes are equal: 

Theorem 4.1 (Persistence equivalence theorem, [10J). Consider per- 
sistent homology of two filtered complexes X and Y . Let (pi : — y 



H*(X ) > H*{Xi) > ... y H*(X n -i) > H*(X n ) 



H*{Y ) > y ... y H*(Y n -i) y H*{Y n ) 



If the (pi are isomorphisms and all the squares commute, then the 
persistence diagrams of X and Y are the same. 

4.5. Computing persistence. Let us have a filtered chain complex 
K.. Boundary matrix d of K encodes the boundary relations between 
cells of different dimensions. Column % corresponds to the boundary of 
cell Cj, row j corresponds to the coboundary of cell Cj. In case of Z 2 
coefficients it can be defined as follows: 



By K,(ci,Cj) := d(i,j) we denote the incidence index of cells c, and 
Cj. In order to compute persistence, a sorted boundary matrix is re- 
quired: For two columns (or rows) i < j, the corresponding cells must 
satisfy: g(ci) < g(cj). Using such a matrix, we can compute persis- 
tence using matrix-reduction Algorithm [l] as defined in [10]. The value 
low(i) marks the maximum (lowest) position of a one in column i, if 
any. We assume it to be zero for zeroed columns. We say that there 
is a collision at column j, if there exists a column k < j such that 



H.(Yi): 




1 if Cj is a face of Cj 
otherwise 



8 



PAWEL DLOTKO AND HUBERT WAGNER 



low(k) = low(j), provided low(k) and low(j) are nonzero. The matrix 
is said to be reduced if there are no collisions. In such all the 

lowest ones are unique. As proven in [10], persistent homology is fully 
determined by the positions of lowest ones in the reduced sorted matrix: 
If column i is zero, corresponding p-dimensional cell q creates a infinite 
p-dimensional persistent homology class. For a non-zero column j with 
k = low(j), the corresponding (p+l)-cell Cj kills a persistent homology 
class created by p-cell cp.- 

The algorithm proceeds with columns from left to right, removing 
any collisions. Later we will analyze the behavior of the algorithm to 
prove the correctness of our simplification algorithm. 



Algorithm 1 Compute reduced matrix 

Input: Sorted binary matrix d of size n x n 
Output: Reduced binary matrix R, which encodes persistence 
1: R:=d 

2: for j := 1 to n do 

3: while there exists in R a nonzero column k < j with low(k) = low(j) 
do 

4: add column k to column j (mod 2) and store as column j 



A simple illustration of Algorithm [T] can be found in Figure [8] in the 
Appendix. 

4.6. Algorithms and their complexity. Applying the presented matrix- 
reduction algorithm to the input complex is the standard way to com- 
pute persistent homology groups. It works for general complexes in 
arbitrary dimensions. The worst-case complexity is 0(n 3 ), where n is 
the size of the input complex. Milosavljevic et al. [28] showed that 
persistent homology can be computed in matrix multiplication time 
0(n") where the currently best estimation of u is 2.3727. Chen and 
Kerber [I] proposed a randomized algorithm to compute only pairs 
with persistence above a chosen threshold. Despite improving the the- 
oretical complexity, it is unclear whether these methods are better in 
practice. 

When focusing on 0-dimensional homology, union-find data struc- 
tures can be used to compute persistence in time 0{na{n)) [TO], where 
a is the inverse of the Ackermann functions and n the size of the input. 

A recent variation of the standard algorithm, introduced by Chen 
and Kerber [5] significantly reduces the amount of computations. This 
idea was also used in [37] to compute persistence for n— dimensional 
images. In general, the regular structure of cubical complexes can be 
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exploited, which allows for handling large inputs. In such a situation 
the size of the boundary matrix is the main obstacle. Preprocessing the 
input complex using discrete Morse theory, as proposed by Robins [M] , 
significantly reduces the size of the boundary matrix, while preserving 
persistnce. In case of 3D grayscale images, an efficient parallel imple- 
mentation was proposed in [16J, allowing for handling large (~ 1200 3 ) 
images on commodity hardware. The standard matrix-reduction algo- 
rithm is used in the final step of computations. 

The approach by Robins works for arbitrary complexes and in di- 
mension three the preprocessing results in the smallest possible bound- 
ary matrix [3lj (counting the number of rows/ columns). The algorithm 
used in [31] depends crucially on simple-homotopy theory, which makes 
it hard to directly generalize the optimality result to higher dimensions. 
A recent paper by Mischaikow et al. [29J proposes a handy theoretical 
framework, where discrete Morse theory is extended from complexes to 
filtrations. 

In our approach, we use the existing algorithms for discrete Morse 
complex construction. We use them to iteratively simplify the input 
complex. Note that our approach is significantly different from the 
simplification scheme by Pascucci et al., where Morse complexes are 
iteratively simplified in terms of (roughly speaking) topology or in a 
sense: persistent homology. Our aim is different: persistence in never 
affected, and the simplification is only in terms of the number of cells 
representing the complex. 



5. Discrete Morse Theory. 

5.1. Morse matching and Morse graph. In this section an algo- 
rithmic introduction to discrete Morse theory is given. For further 
theoretical details please consult [131 [25] . Let us have a complex /C. 
Discrete Morse theory partitions the cells of K, into matched cells and 
critical cells. The critical cells, together with a boundary operator we 
describe later, form a chain complex called the Morse complex. Impor- 
tantly, this Morse complex has homology isomorphic with the homology 
of the initial complex. A procedure to compute Morse complex is given 
in Algorithm [2j 

The first step in constructing a Morse complex requires finding an 
acyclic Morse matching M of the complex /C. The matching is a partial 
map M : K. -+> K.. Each cell of K. can be matched with exactly one of 
its co-faces. Some cells can remain unmatched, and are called critical. 
These cells constitute the resulting Morse complex. 
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Algorithm 2 Compute Morse complex 

Input: Input complex tC 
Output: Resulting Morse complex 

1: M := acyclic Morse matching on K, 

2: C := list of critical cells of M 

3: G := Morse graph of M, C 

4: bd := compute Morse Boundary of G (as in Algorithm [4] 
5: return (C, bd) 



Let us introduce a concept of a Morse graph of a complex K, and 
matching M. It is a directed graph whose vertices are formed by cells 
of a complex. A directed edge from vertex A to vertex b is added 
whenever b is in the boundary of A and the cells A and b are not 
matched in M. If they are matched in M, a direct edge from b to 
A is added in the Morse graph. The matching M is called acyclic 
if the corresponding Morse graph is a directed acyclic graph (DAG). 
The paths of this graph are often called V— paths. There are various 
strategies of obtaining acyclic Morse matchings. Later in this paper we 
assume that every matching is acyclic. 

A Morse matching is called perfect if each unmatched cell corre- 
sponds to a homology generator of the original complex. (We choose 
to talk about (perfect) matchings, but in the literature (perfect) Morse 
functions, vector-fields and complexes are discussed.) This depends on 
the choice of coefficients, for example in case of a field coefficients F, 
we can talk about F-perfect matchings Some spaces do not admit 
perfect matchings, for example the Dunce hat (presented in Section [I]), 
being contractible but non-collapsible. In this case some critical cells 
are spurious, in the sense that they do not correspond to any homology 
generators. 

One can try to construct a best possible matching, minimizing the 
number of critical cells. This problem is known to be NP-complete 
and MAXSNP-hard [22] • It means that computing the best possi- 
ble matching is computationally expensive and there is no hope to 
find a fully-polynomial time approximation strategy. As a result, no 
polynomial-time algorithm can give arbitrarily good bounds on the 
number of spurious critical cells. 

Once an acyclic Morse matching M is obtained, we proceed with 
computing the Morse boundary. This procedure is described in [13]. 
The idea is illustrated in Figure [3] Forman [13] proved that the re- 
sulting Morse complex has isomorphic homology to the homology of 
the initial complex. He also provided a formula which computes the 
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boundary of each cell in a Morse complex. Kozlov generalized these 
proofs to the setting of arbitrary chain complexes [25]. One can use 
Morse complex construction to compute homology. Later we show that 
similar construction can be also used to compute persistence. 

5.2. Discrete Morse theory for filtered complexes. Recently there 
were first successful attempts to use discrete Morse theory to compute 
persistence [HU [TO] (in case of 3-d gray-scale images), [2] (in case of 2- 
manifolds). The first successful attempt to provide a Morse-theoretic 
categorical framework for persistent homology was made in [29J. 

In this section we will first recall the basic ideas from [29]. We say 
that the Morse matching M is compatible with filtration of K. if for every 
matched A 6 fC, g{A) = g(M(A)]^ In other words, the matchings 
are made between elements of the same filtration level. Consequently, 
directed paths cannot move upwards the filtration (if they did, we 
would lose the sub-complex filtration property in the corresponding 
Morse complex, because a cell could enter the filtration strictly before 
its faces). 

The key result in [29J is that the persistence diagram of a filtered 
complex /C and a Morse complex M(/C) with the Morse matchings 
compatible with filtration are the same. In Figure [2] an example of 
Morse matchings compatible and non-compatible with filtration are 
shown. 




1 



Figure 2. On the left the Morse matching compatible 
with filtration. In this case persistent homology for both 
initial and the Morse complex is [0, oo] in dimension 
and [3, oo] in dimension 1. On the right a correct Morse 
matching in a sense of standard Discrete Morse Theory 
which is not compatible with filtration is depicted. The 
persistent homology of the Morse complex on the right 
in dimension one is [1, oo] which is not correct. 



2 By M(A) we denote the element matched with A. 
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5.3. Computing Morse complex with graph algorithms. In this 
section we will show that the entire Morse complex construction can 
be computed using standard graph algorithms. The chain complex 
(with Z2 coefficients) can be interpreted as a graph - namely the Morse 
graph. We assume that initially no matchings are made. The whole 
construction can be divided into two essential parts: finding an acyclic 
Morse matching and computing the Morse boundary. Both parts are 
described below. We want to point out that for presentation's sake the 
algorithms presented in this section are not necessarily optimal. For 
that reason we also present just a version for binary coefficients. They 
can be easily generalized to arbitrary field coefficients. 




Figure 3. This picture shows the process of comput- 
ing the Morse boundary of cell A. A Morse matching 
is shown. The part of corresponding (acyclic) Morse 
graph is build. Finally boundary relations between cells 
are computed, by following paths emanating from the 
boundary of cell A and ending at critical cells. Note 
that there are two paths between cells A and d, which 
yields boundary relation equal to zero for coefficients in 
Z 2 (d is not in the boundary of A). 
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5.3.1. Computing acyclic matching. There are many possible graph- 
based strategies for Morse matchings. Some strategies based on the 
idea of BFS spanning tree have been described in the literature j22j 
|26| [38] • In general, the problem of performing a Morse matching is 
equivalent to the following one - which directed edges in Morse graph 
can be reversed so that the graph remains acyclic. This is closely 
related to the minimum feedback arc set problem. While this problem 
is NP-hard, efficient approximation schemes exist [T2] . 

For an illustration, a basic algorithm descried in [T5] will be reminded 
in Algorithm [3} The proof that the obtained graph is acyclic is easy 



Algorithm 3 Compute Morse matching 
Input: Morse graph G; 

Output: Changed graph G (edges between matched elements are reversed); 

1: while Not all vertices of G are marked do 

2: Let i be the minimal dimension of unmarked element in G\ 

3: Pick A, unmarked i-dimensional cell in G and mark as critical] 

4: while There exists unmatched element A € G with unique unmatched 

element b in boundary do 
5: Make a Morse matching (A, b), i.e. reverse an edge from A to b in 
G- 

6: Mark A and b as matched; 



but technical. Therefore it will not be presented here. 

To put restrictions on the matching one should modify line [4] of 
Algorithm |3j In particular, one can ensure that the matchings are 
compatible with filtration. This is needed in case of persistence com- 
putations. 

5.3.2. Computing Morse boundaries. Let us state the problem of com- 
puting Morse boundary in terms of graph theory. The Morse boundary 
of a critical cell is formed by the set of critical cells which are reachable 
by an odd number of paths in the Morse graph. This is a special case 
of Forman's formula [13J in case of Z 2 coefficients. Since the number of 
such paths can grow exponentially in the number of cells, brute force 
calculation is ineffective (as noted in [29]). However, this problem can 
be solved efficiently, exploiting the fact that the Morse graph is acyclic. 

The following Algorithm [4] is simple to implement, and works in 
pessimistic time 0(c* (\\V\\ + \\E\\)), where V and E are respectively 
the vertices and edges of graph G and c is the number of critical cells. 
For a runtime- and memory-optimal one see [16J. 
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Let P s {t) denote the number of distinct paths leading from vertex s 
to t, and prev(v) = {x \ (x,v) G E} is the set of vertices preceding v 
in the directed graph. We have an obvious recurrence relation: 



u 



1 for u = s 

J2veprev(u) P s( V ) for U ^ S 

This recurrence can be computed directly and also efficiently using 
memoization, but we propose an elegant graph-theoretical algorithm. 



To compute P s (v) the summation is done indirectly in line 10 of the 
Algorithm [4] by adding the value P s (u) where u E prev(v). 



Algorithm 4 Compute Morse boundaries from Morse graph 

Input: Directed Morse graph G := (V,E) 
Output: Boundary relation d of the Morse complex 
1 

2 
3 
4 
5 
G 
7 



9 

10 



sort G topologically 
for each critical vertex s do 

assign P s (v) := for each vertex v 7^ s 

assign P s (s) := 1 

for each vertex c following s in topological order do 
if c is critical then 

d(s,c) := P s (c) mod 2 
else 

for each v \ (c, v) E E do 

P s {v) += P s (c) 



Theorem 5.1. Algorithm^is correct. 

Proof. We prove the correctness of the algorithm by induction on the 
iteration of the loop in line [5] The desired invariant is that whenever 
P s {c) is used, this value is already final and correct. Clearly, the value 
P s (s) is initialized correctly in line|4| so the correct value is used during 
the first iteration. Assume that the invariant holds for the first % itera- 
tions and vertex c is now processed. Note that the value of c depends 
on the values of prev(c). Since we proceed in topological-sort order, 
all the vertices in prev(c) have already been processed. By inductive 
assumption the values used to compute P s {c) were correct and final, 
therefore the value P s {c) is also correct and final. □ 

6. Iterated Morse Complex for homology 

In this section the concept of iterated Morse complex is presented. 
Normally one aims at finding a Morse complex minimizing the number 
of critical cells. As mentioned earlier this is a hard algorithmic problem 
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and we do not tackle it. Instead we use an algorithm to iteratively 
construct a sequence of Morse complexes. If at a certain stage the 
obtained Morse complex is far from optimal, further iterations will be 
necessary to compute the homology of the considered complex. Still, 
the worst case computational time is cubical. The results presented in 
this section has already been sketched in [38J. 

Let C be a category of chain complexes and let M : C — > C be 
a functor taking a chain complex and assigning it a Morse complex 
constructed on it. There are many possible strategies to construct 
Morse complexes. We assume that if a chain complex /C G C has 
some available Morse matchings, M does at least one of them^j This 
property of M will be referred to as vitality. Except from vitality no 
extra assumptions are put on M. 

For a given chain complex K. E C, iterated Morse complex M°°(/C) is 
the fixed point of the iteration M(/C),M 2 (/C) = M(M(/C)), M 3 (/C), . . .. 
It is clear that ||/C|| > ||M(/C)|| > ||M 2 (/C)|| > . . .. Moreover due 
to the vitality of M the above inequalities are strict as long as there 
are some Morse matchings to be made in the intermediate complexes. 
Therefore, the fixed point M°°(/C) is obtained in a finite number of 
iterations. Below we show that M°°(/C) gives an instant information 
about homology of /C. To achieve this, we will use algebraic version of 
discrete Morse theory due to Kozlov [25]. It states that two elements 
A, B can be matched if and only if k(A, B) is invertible. Since in 
this paper we consider only homology with field coefficients, k(A, B) is 
always invertible provided it is nonzero. This fact implies the following 
straightforward lemmas: 

Lemma 6.1. For every A e M°°(/C) both boundary and coboundary of 
A are empty. 

Lemma 6.2. = \\{A e M°°(/C)| dim A = 

The proof of the first lemma is a direct consequence of vitality of 
M. If there exists A 6 M°°(/C) with B in (co)boundary, then M would 
eventually make a Morse matching between A and B. 

The second lemma is a direct consequence of the first one. Once every 
element has empty boundary, it is a cycle. Once it has empty cobound- 
ary, it cannot be a boundary. Therefore every element in M°° (K) gen- 
erates a homology class of /C. Moreover, due to Section 11.3 in [25] it is 



The simplest example of algorithm that fulfill this requirement is M which 
searches for a first possible Morse matching in /C, makes it, and computes the 
Morse complex. 
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clear that the homology of MP(/C) and MP +1 (/C) are isomorphic. There- 
fore the homologies of the complex are preserved through the entire 
iteration. 

As already mentioned, the idea of iteration of Morse complex con- 
struction implies that we do not have to construct near optimal Morse 
complexes. Let n be the cardinality of /C. In the worst case, after \n/2\ 
iterations of Morse complex procedure, an iterated Morse complex is 
obtainecFJ This is a consequence of vitality of M. 



Algorithm 5 Compute homology with iterated Morse decomposition 

Input: Initial complex C of dimension d 
Output: Betti numbers 
1: while true do 

2: M := Build Morse complex of C (Algorithm [2]) 
3: if M = C then 
4: break 
C := M 
5: for i := to d do 

6: := number of i-dimensional cells of C 



Algorithm [5] describes how to comute the Betti numbers using it- 
erated Morse decomposition. An example of the Morse complex con- 
struction on a Dunce hat has already been presented in Section [T] In 
case of the Dunce hat there does not exist a perfect Morse complex. At 
the end of this section, in Figure [5j we show a simple example of the 
presented construction, using a sub-optimal algorithm M to construct 
Morse complexes with sub-optimal Morse matchings. 

7. Iterated Morse Complex for persistent homology 

In [29] it is shown that when a Morse complex is constructed based 
on a Morse matching compatible with filtration, persistent homology of 
the initial complex and the one of the Morse complex are isomorphic. 
Here we provide a further consequence of this result. Let us take a 
vital functor M acting from a category of filtered chain complexes to 
itself. We assume that the Morse matching used to construct M(/C) is 
compatible with filtration of /C. Filtration values of cells in M(/C) are 
inherited from the filtration of cells in K,. As in Section [6j we construct 
M°°(/C). 



Maximal number of iterations needed is the floor of cardinality of basis of /C 
divided by two minus sum of all the Betti numbers of K,. 
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In this section we show that each cell in M°°(/C) either creates or 
kills a feature of nonzero persistence. Therefore, if we want to mini- 
mize the number of cells, the resulting complex is the minimal complex 
encoding persistence of the original complex. An example is presented 
in Figure |4j 




Figure 4. (a) The initial complex K, with a few Morse 
matchings indicated with arrows, (b) Morse complex 
M 1 (/C) obtained from /C. A few possible Morse matchings 
indicated with arrows, (c) Final Morse complex M 2 (/C) 
obtained from M 1 (/C). 

In it is only assumed that the filtered chain complex is given at 
the input. Since at each stage of iterated Morse complex computation 
we have a chain complex, one can iterate further the construction. 
The main strength of our approach is based on the following novel 
observation. The resulting complex M°° (/C) has the following property: 
For every A G M°° (/C) and for every b±, . . . , b n in boundary of A we have 
g(A) > g(bi), . . . , g(b n ). It is because if there existed bi in the boundary 
of A such that g(bj) = g{A), then a Morse matching could be made 
between A and 6, (since coefficients in a field are used). This would 
contradict the vitality assumption of M. Having this simple property 
of the resulting complex M°° (/C) we can now present the main theorem 
of this section: 

Theorem 7.1. Let M°°(/C) be the iterated Morse complex obtained 
from the initial filtered chain complex K. by iterative construction of 
Morse complexes using Morse matchings compatible with filtration. Then 
K. and M°°(/C) have the same persistence. (Correctness) 
Every element A 6 M°°(/C) either starts or terminates a nonzero length 
persistence interval. (Optimality) 

Proof. Correctness: To show that the algorithm is correct we will apply 
the Persistence Equivalence Theorem for each iteration: 
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bd M 



... ^(Mj +1 ) #*(Mj+i) ... 

We need to show two things: 

(1) That the vertical maps are isomorphisms on homology level. 

(2) That all squares commute. 

First note that the vertical maps send each chain in the input com- 
plex to a corresponding chain of the Morse complex. This is equivalent 
to computing Morse boundaries, as described in Section |5j We remind 
that to find a corresponding Morse chain we follow appropriate V-paths 
in the Morse graph. 

1) Vertical arrows are isomorphisms: this is a consequence of The- 
orem 2.1 from [25], which states that the upper chain complex M is 
decomposed by the Morse construction into an acyclic part and the 
Morse complex having homology isomorphic with M. 

2) To prove that the squares commute for each i, I, let us take a chain 
c g Mj = cj. We show that (ftdftj o j z )(c) = (Jfc, o 6df)(c). 

Down and right {ki o bdf 1 ): If Cj is critical, it is unchanged by the 
vertical map. Otherwise, we follow the paths of the Morse graph to 
compute the corresponding chain in the Morse complex. Repeating this 
computation for every Cj, the value of baj 1 on c is obtained. Moving 
right with inclusion, the chain remains the same. 

Right and down (bd^ 1 o first the chain c is inserted by inclusion 
into level I + 1 of filtration, so it is unchanged. But now we move with 
the vertical arrow, which might be richer on this level, as additional 
paths enter the Morse graph. Note that since we force the paths to be 
non-increasing with filtration, bdf^_ 1 restricted to level I is the same as 
bd^ 1 . In other words, any V-path starting at Cj at level at most I can 
only reach cells of lower or equal filtration values. In particular, it will 
never reach any critical cell introduced at level l + l. 

Therefore the two images of chain c are the same and the diagram 
commutes. 

We can apply Persistence Equivalence Theorem and finish the proof 
that persistent homology is unchanged during our iterated Morse com- 
plex construction. 

Optimality: 

We want to show that the simplified complex, M°°(/C) contains only 
significant information, that is no intervals of persistence zero are 
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present. The argument is based on the analysis of the behavior of the 
standard (left-to-right) matrix-reduction algorithm (Algorithm [T]) run 
on the boundary matrix of the final iterated Morse complex M°°(/C). 
The lowest-ones of the reduced matrix directly indicate persistence in- 
tervals. Zero columns indicate that a given cell creates an infinite 
interval. 

The argument is inductive with respect to the iteration of the outer 
loop in the Algorithm [T] Specifically, we consider the first k reduced 
columns of the matrix. For k = this submatrix is empty, so there are 
no zero-persistence pairs. 

Let us assume that the argument holds for some k > 0. 

Suppose by contradiction that there is a zero-length persistence in- 
terval generated by a pair (A, b), where A is the k + 1 column. It means 
that at this stage we have the following situation (dots mark arbitrary 
entries): 

The matrix after k + 1 iterations (first k + 1 columns are reduced). 

k + 1 
A 

b 0... 1 




The matrix after k iterations (first k columns are reduced). 

k + 1 
Al AO 



b 



61 1 1 

There are two possibilities: 

(1) During the reduction of AO there were no collisions, so ^40 = 
A is a cell in the complex M°°(/C). Then, since g{A) = g{b) 
and k(A, b) ^ 0, it was possible to make a Morse matching 
between A and b, which gives a contradiction with the fact that 
all possible Morse matchings compatible with filtration were 
made in M°°(/C). 
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(2) A is represented as a sum of the preceding columns, which gen- 
erated collisions. In this case A = AO + Al + . . . An, for AO 
being the column in the unreduced matrix and Al, . . . , An be- 
ing columns preceding AO in the matrix. 

For the proof we need to find the lowest nonzero position of all 
the columns A4, i G {1, . . . ,n}. During the process of reducing 
a single column the lowest-ones can never increase, therefore 
the first collision yields the lowest-one we search for. We call 
the column Al and this lowest position bl and note that g(b) < 
g(bl). Also note that bl marks the lowest one in the reduced 
column Al and the original column AO, so g(bl) < g{Al) and 
g(bl) < g(A0) = g(A). 

From the assumption we have g(A) = g(b). From the filtra- 
tion of the complex we have g(Al) < g(A) and g(b) < g(bl). 
Putting this together we get: g(A) = g(b) < g(bl) < g(Al) < 
g{A), consequently g(bl) = g(Al). This means that there was 
a zero-persistence pair within the first k columns. This contra- 
dicts the inductive assumption. 

□ 



We have the following theorem which is a direct consequence of The- 
orem I7.lt 

Theorem 7.2. Let /C be the initial filtered chain complex. Let p be the 
number of finite and k the number of infinite persistence intervals of 
K. Then ||M°°(/C)|| =2p + k. 

This theorem indicates that the complex M°° {JC) is the minimal com- 
plex encoding the persistence of the initial complex. We can now ap- 
ply the matrix reduction method, to get persistence intervals in time 
0((2p + k) 3 ). Therefore if the number of persistence intervals is small, 
this computation can be efficient. As an alternative, we propose a new 
algorithm presented in Section [8j which relies only on graph operations. 

Algorithm [6] describes the simplification procedure. To compute per- 
sistence based on simplified complex C use Algorithm [T} We want 
to point out that there is no obvious way of relaxing the condition 
g{A) = g(M(A)) for matched elements. See the Appendix for more 
details. 

At the end of this section, in Figure [5] we present an example of the 
iterated Morse complex construction. 
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Algorithm 6 Simplification for persistence computations 

Input: Initial filtered complex C of dimension d 
Output: Simplified complex C, having the same persistent homology 
1: while true do 

2: M := Build Morse complex of C using only matchings compatible with 

filtration 
3: if M = C then 
4: break 

C := M 



8. Persistence intervals via iterated Morse approach. 

In this section we compute persistence intervals using the iterated 
Morse complex approach. It will be shown that the presented approach 
can be interpreted as a variation of matrix-reduction algorithm. It is 
however based on graph theory, rather than matrix algebra. 

We stress that in this section, unlike previous one, we allow to make 
Morse matchings between elements having different level of filtration. 
Therefore, we will construct Morse matchings that are not compatible 
with filtration in a sense that element A can be matched with b if 
9(A) > g(b). 

Algorithm [7] is used to compute persistence intervals of M°° (JC) . 

We want to point out that the Morse boundary procedure is always 
performed on the whole complex M°°(/C) even if the matchings are 
made on a proper subcomplex M°°(/C). We start from the complex 
Mj D (/C). Since M°°(/C) is an iterated Morse complex, is clear that 
in Mf^K.) and Mf(K) \ M^fX) there are no Morse matchings to 
be made. But in Mj° (/C) there can be a Morse matching (A, b) such 
that A G Mf(K.) \ M^/C) and b G M^fX). This means that b G 
Mf^K.) is a homology generator in Mf^K.) which is killed in Mf(JC) 
(since the matching (A, b) can be made without changing homology 
of Mj D (/C)). Making such a matching indicates a persistence interval 
generated by the pair (A, b), since the homology class generated by b 
is killed by A. We assume that all possible matchings in My°(/C) are 
made (by iterating Morse complex construction) before proceeding to 
Mf +1 ()C). 

When processing complex M?°(/C) we assume that in M°^ 1 (/C) no 
more Morse matchings can be made. We are searching for matchings 
(A, b) such that A G Mf (K) \ M°^(/C) and b G M^/C). If two or 
more elements b%, . . . ,b n can be matched with A we always choose be- 
having maximal value of filtration. Morse matching (A, b) indicates an 
interval generated by (A, b), since b generates nontrivial homology class 
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Figure 5. On the top left the initial filtered complex K, 
is depicted. On the top right, the first iteration, M 1 (/C) of 
the Morse complex compatible with filtration construc- 
tion. On the middle left with red arrows the second iter- 
ation M 2 (/C), and on the middle right the third and final 
iteration of M°°(/C) construction. On the bottom left, 
the cells of M°° (/C) are named, and on bottom right the 
boundary relation is presented in a form of diagram. The 
levels indicate the gradation - vertices a, b at the bot- 
tom, edges c, d, e in the middle and faces X and Y at the 
top. The numbers indicate filtration values of elements 
and arrows - boundary relation. 



in the level of g(b), which becomes trivial or identical to some other 
class born earlier in the level of g(A). 
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Algorithm 7 Compute persistent homology via Iterated Morse com- 
plex. 

Input: Filtered iterated Morse complex C = M°°(/C). 
Output: Persistence intervals of M°°(/C). 

1: S := empty multiset of persistence intervals; 

2: int / := second filtration level of C; 

3: int / := last filtration level of C; 

4: M is a vital strategy to make Morse matchings. Matchings between 
elements of different nitrations are allowed in M. If element A can be 
matched with two (or more) elements b%, b2---, we match it with the one 
with maximal filtration value. 

5: for int i := f to I do 

6: while true do 

7: Construct Morse matching on Ci using M (remark: for every (A, b) 

matched by M we have g(b) < g(A) = i.). 
8: if Nothing was matched by M then 
9: break 

10: for Every (A, b) matched by M do 
11: S:=SU[g(b),g(A)]; 

12: C := Morse complex on C constructed based on Morse matching 
M; 

13: for every cell c in complex C do 
14: S := SU[g(c),oo\; 
15: return S; 



When all the possible matchings are made, in the last for loop the 
unmatched elements are found. They generate infinite persistence in- 
tervals. 

It is clear that the levels of filtration need to be processed in order. 
See the Appendix for more details. 

Let us now show that the presented technique can be interpreted in 
therms of the standard algebraic Algorithm [TJ 

Theorem 8.1. Let K, be a filtered chain complex. Let M°°(/C) be the 
iterated Morse complex described in Section [6| Then the persistence 
intervals of K, and the intervals obtained from M°°(/C) by the described 
algorithm are the same. 



Proof. From Theorem 7^ it is clear that the persistence intervals of 
K. and M°°(/C) are the same. Let (A,b) be the first Morse matching 
made by the presented algorithm (i.e. there does not exist a possible 
matching (A', V) such that g(A') < g(A) and there does not exist b 1 , a 
face of A with g(b') > g(b)). 
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Let M (Aib) (M 00 (/C)) be the Morse complex obtained from M°°(/C) by 
making a Morse matching (A, b). To prove the theorem it suffices to 
show that the multiset of persistence intervals of M°°(/C) minus the 
interval [g(A), g(b)} is equal to the multiset of persistence intervals of 
M (A6) (M°°(/C)f] ' 

First let us remind how the procedure to compute Morse boundary 
works and how it is interpreted in matrix-reduction algorithm. The 
details can be found in [13J. Let us assume just one Morse matching, 
(A, b), is made. And also that b is in the boundary of A,Ai,..., A n . 
Then boundaries of Ai,...,A n need to be changed in the following 
way: dAi = dAi \ b U dA \ & i.e. b is replaced in boundary of A% with 
boundary of A excluding b, as in Figure [6] 




Figure 6. Illustration how performing a single pairing 
changes the complex. 

Now, suppose Algorithm [I] is run on the complex M°°(X). Without 
loss of generality we may assume that the column A is the first nonzero 
column in the matrix. Since there cannot be a collision there, Algo- 
rithm [l] leaves the column A unchanged and the interval [g(b), g(A)] is 
obtained in dimension of b - as in the case of the algorithm presented 
in this section. But, row b may cause some collisions later in the course 
of execution of Algorithm [T] Suppose the first collision in Algorithm [T] 
occurs: 

A A, 

b ... 1 1 




To remove this collision, Algorithm [T] performs the addition Ai : = 
Ai + A. In the algorithm presented in this section, when processing 

5 This follows from the fact that Mau in , ,(A„,6„) = Mu^u o ... o y^(A n ,b n ) ■ 
Easy proof is left for the reader. 
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cell A the matching (A, b) was made and the Morse boundaries were 
computed. As one can see, the computations of Morse boundary after 
the matching (A, b) is simply equivalent to summing Ai := Ai + A for 
all Ai having b in boundary. Therefore all the future collisions caused 
by b are resolved. Consequently making a Morse matching is equivalent 
to resolving all the future collisions at once. This simple observation 
proves the theorem. □ 

In Figure [7] an illustration of the presented procedure on the complex 
M°°(/C) from Figure [5] is given. 

We want to point out that this approach promises to parallelize well. 
The details will be presented in a more technical paper [U] . 

Moreover, the approach will be as scalable as the implementations 
of the underlying graph algorithms. Therefore, using the available, 
mature libraries, we hope to achieve good practical performance. 

We also want to point out that Algorithm [7] can be used with minor 
modification for initial filtered complex /C. Easy changes are left for 
the reader. 

9. Conclusions 

In our opinion the presented technique has several advantages, com- 
paring to the standard way of computing homology and persistence: 

(1) It is combinatorial, does not require any algebraic matrix oper- 
ation. 

(2) It is based on graph theory - we can use efficient algorithm 
(exact and approximate) and their existing implementations - 
in particular libraries for distributed graph operations [27J. 

(3) It is intuitive - it is easy to visualize the process of homology 
computations. 

We are aware of some drawbacks and complications of our method: 

(1) There may exist bad cases, where the complexity will be unsat- 
isfactory. 

(2) This technique might not be suitable for cubical data. Existing 
methods [161 EZ] rer y on the compact representation of cubical 
grids. In our case the complex need not be a cubical complex 
after the first iteration, preventing us from storing it efficiently. 

The presented techniques can be used to formalize and generalize 
homology-preserving properties of graph pyramids used in image recog- 
nition [31] . Moreover they can be beneficial in verified homology com- 
putation |20j by avoiding matrix operations which are costly to verify 
automatically. 
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Initial complex M_ 

(Y,6) (X,6) 



g -bi (c,4) (d,5) (e,3) 



(a,1) (b,1) 



Complex after second iteration 

f(Y6K (X,6) 




(a,1) 



Complex after first iteration 

(Y.6) (X6) 

Y 

[a d J (c^4) (d*5) Je, 

x I \£ U 

.^^^J (a,1) (b,1.) 



Complex after third iteration 




(a,1) 



Figure 7. First complex on the top left - complex 
M°°(/C) from Figure [5j The complex on the top right 
- on filtration value 3 the first possible Morse matching 
between b and e appears (indicated on the left with ar- 
row, and with ellipse on the right) is made. When the 
matching is constructed, the persistence interval [1, 3] in 
dimension zero is reported. Third complex on the bot- 
tom left is obtained as a result. There, on the level 6 
a matching between Y and d can be made (we want to 
point out that there is also possible matching between 
Y and c, however filtration value of d is higher than fil- 
tration value of c). After the matching, the persistent 
interval [5, 6] in dimension one is reported. The final 
complex on the bottom right- the last possible match- 
ing between X and c is made. After the matching, the 
persistent interval [4, 6] in dimension 1 is reported. The 
remaining cell in the complex is the vertex a it corre- 
spond to infinite persistence interval [1, oo] in dimension 
zero. 



Summarizing, in this paper a novel technique of homology and per- 
sistent homology computations has been presented. It indicates that 
problem of computing (persistent) homology can be solved by itera- 
tively applying standard graph algorithms. We hope that this approach 
will lead to a scalable implementation for (persistent) homology com- 
putations. 
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10. Appendix 

10.1. Example of matrix reduction computations. A simple ex- 
ample of persistence intervals computations with the Algorithm [T] are 
presented in Figure |HJ On the left, the initial complex. We assume 
that the filtration value for every vertex is 0. The filtration value of 
edges are given in the picture. On the upper left, the initial boundary 
matrix. As one can see, the only collision is between columns cd and 
bd. Therefore we have column cd = cd + bd on the upper right. Then 
a collision between cd and ac appears and we set cd = cd + ac on the 
lower left. There again we have a collision cd with ab which is removed 
in the lower right by setting cd = cd + ab. On the matrix in lower 
right there are no more collisions, therefore we can read persistence 
intervals out of it. Lowest one in column ac indicates that edge ac kills 
connected component created by c, which gives an interval [0, 1] in di- 
mension 0. Analogously lowest one in column bd induces an interval 
[0, 1] in dimension 0. Lowest one in column ab indicates that edge ab 
kills a connected component created in b, which gives an interval [0, 2] 
in dimension 0. Zero column cd induces an infinite interval [3, oo] in 
dimension one. 

10.2. Relaxing the tolerance. In Figure [9] we show that if one makes 
a Morse matchings between elements (A,M(A)) such that \g(A) — 
g(M(A))\ < e, one may get arbitrarily large differences in the output 
persistence intervals. 



10.3. Processing order. In Figure 10 it is shown what happens if we 



do not proceed in Algorithm [7] in the filtration order. 
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